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Abstract. For stochastic implicit Taylor methods that use an iterative scheme to com- 
pute their numerical solution, stochastic B-series and corresponding growth functions are 
constructed. From these, convergence results based on the order of the underlying Taylor 
method, the choice of the iteration method, the predictor and the number of iterations, for 
Ito and Stratonovich SDEs, and for weak as well as strong convergence are derived. As 
special case, also the application of Taylor methods to ODEs is considered. The theory is 
supported by numerical experiments. 

1. Introduction 

Besides stochastic Runge-Kutta methods, one important class of schemes to approxi- 
mate the solution of stochastic differential equations (SDEs) are stochastic Taylor methods. 
As in the deterministic setting [1], they are especially suitable for problems with not too 
high dimension, and here especially in the case of strong approximation, because weak 
approximation of low-dimensional problems can often be done more efficiently by nu- 
merically solving the corresponding deterministic PDE problem obtained by applying the 
Feynman-Kac formula. For solving stiff SDEs, implicit methods have to be considered, as 
is illustrated in the following two examples. 

Example 1.1 (see [12]). Consider the linear Ito-SDE 

(1) dX{t)=^X{t)dt + aX{t)dW{t), X{0)=xo, 

with |i , C7 G C. We assume that the exact solution is mean-square stable, i. e. 

limE(|X(f)|^) ==0, 

which is the case if and only if l^ifj. + |(7p < 0. To achieve that also the numerical ap- 
proximation Y„ obtained with the (explicit) Euler-Maruyama scheme with step size h is 
mean-square stable, i.e. lim„^<joE(|y„|^) = 0, we have to restrict the step size accord- 
ing to h < ho := — ^^I'^j"^ , whereas for h > ho the numerical approximations explode, 

lim„_).c„E(|F„p) = °°. In contrast to this, the semi-implicit Euler scheme is mean-square 
stable without any step size restriction. For a numerical affirmation, see Figure 1. 

Example 1 .2. Consider the following stochastic Van der Pol equation, 
dXi{t) ^X2{t) dt, 

dX2{t) = (^(1 -Xi{tf)X2{t) -Xi{t)) dt + e{i ~Xi{tf)X2{t) dw{t), 

Xi (0)^X0,1, X2(0)=xo,2. 
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(a) Euler-Maruyama scheme (explicit) (b) Semi-implicit Euler scheme 



Figure 1. Approximation results for the linear test equation (1) with 
= —3, a — \/3, and xq = 1 by Euler-Maruyama (explicit) and semi- 
implicit Euler scheme. Here, E(|F(f)p) is approximated as mean over 
10^ simulations. The explicit scheme is only stable for appropriate step 
sizes, the semi-implicit scheme is stable for all step sizes. 




Figure 2. Approximation of Van der Pol equation with /i = 10, = 1, 
xq \ — 2, and jcq 2 = by the explicit and semi-implicit Milstein scheme 
(using the same Brownian path) with step-size h = 0.05. The explicit 
scheme suffers from heavy stability problems and aborts. 

Application of the explicit Milstein scheme (see Example 1.3 with a = and j3 = 0) with 
step-size h — 0.05 to approximate a solution path leads to an explosion of the approxima- 
tion, see Figure 2(a), whereas application of the semi-implicit Milstein scheme, given by 
substituting goiY„) by go(Jn+i) (CC = I and ji = in Example 1.3), yields (for the same 
Brownian path) the result of Figure 2(b). 

Implicit stochastic Taylor methods have been considered both for strong [15, 17] and 
weak [15] approximation. For these methods, the approximation values are only given 
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implicitly. However, in practice these implicit equations are solved by iterative schemes 
like simple iteration or Newton iteration. The "exact numerical" solution can be written 
in terms of B-series [8]. As we will prove in this paper, so can the iterated solution. 
Moreover, for each iteration scheme in question, we will define a growth function. Briefly 
explained, when the exact numerical and the k times iterated solutions are both written in 
terms of B-series, then all terms of these series for which the growth function has a value 
not greater than k coincide. Thus the growth functions give a quite exact description of the 
development of the iterations. B-series and corresponding growth functions for iterated 
solutions have been derived for Runge-Kutta methods apphed to deterministic ordinary 
differential equations [14], differential algebraic equations [13], and SDEs [7]. Somewhat 
surprisingly, the growth functions are exactly the same in all these cases, and, as we will 
show in this paper, this also holds for impUcit Taylor methods. 

The outline of the paper is as follows: First, we will give the SDE to be solved and the 
iterated Taylor methods used for its approximation. In Section 2 stochastic B-series are 
introduced and some useful preliminary results are presented. The main results of the paper 
can be found in Section 3, where the B-series of the iterated solutions are developed and 
the before mentioned growth functions derived. In Section 4, these findings are interpreted 
in terms of the order of the overall scheme, giving concrete results on the order of the 
considered methods depending on the kind and number of iterations, both for SDEs and 
ODEs. Contrary to the results obtained for Runge-Kutta methods [7], the order of the 
iteration error is shown to be independent on whether Ito or Stratonovich SDEs, weak 
or strong convergence is considered. Finally, in Section 5 we present several numerical 
examples to support our theoretical findings. 

Let {fl,s^ ,^^) be a probabihty space. We denote by {X{t))tei the stochastic process 
which is the solution of a J-dimensional SDE defined by 

m 

(2) dX{t)=g(,{X(t))dt + Y.Si{nt))^dWi{t), X{to)=xo, 

1=1 

with an m-dimensional Wiener process {W{t))t>o and/= [to,T]. As usual, (2) is construed 
as abbreviation of 

(3) X{t)=xo+ / 8o{Xis))ds+Y, / gi{X{s))*dWi{s). 

J to J to 

The integral w. 1. 1. the Wiener process has to be interpreted e. g. as Ito integral with -kdWi{s) = 
dWi (s) or as Stratonovich integral with -kdWi (s) = odWi {s) . We assume that the Borel-mea- 
surable coefficients gi -.M.^ ^ are sufficiently differentiable and chosen such that SDE 

(3) has a unique solution. 

To simplify the presentation, we define Wo{s) = s, so that (3) can be written as 

m 

(4) X{t)=xo + '£ / gi{X{s))*dWi{s). 

1=0 -^'0 

Let a discretization = {?o,?i, . . . ,1^} with to < ti < . . . < = T of the time interval / 
with step sizes h„=t„+i— tn for n = 0,1,..., N— I he given. Now, we consider the general 
class of stochastic Taylor methods given by Yq = xq and 

(5) Yn+i =B{^ex,Yn;hn) + B y„+i;/j„) 

forn = 0, 1 , . . . ,Ar - 1 with Y„ = Y{t„), tn e /^ ^ex{^) = 1, 3>,;„(0) = 0. 
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Example 1.3. Consider the family of Milstein schemes applied to an ltd SDE with one- 
dimensional noise, 



(6) Y„+i=Yn + hn{{l-a)go{Yn) + ago{Y„+i)) 

+/(i),ft„((i - j3)^i(y„) + j3gi(y„+i)) + - /3/fi),,„)kUi](i'„). 



and the parameters a,j3 G [0, 1] indicate the degree of implicitness. When a = ^ = Owe 
have the explicit Milstein scheme, with a ^0, ji = a semi-implicit scheme. In all cases, 
the method (6) can be written in the form (5) with 



The terms and <!>,„, refer to the time- and method-dependent part of each term: In this 
case 



The notation will be explained in detail in Section 2. 

What the general method (5) concerns, application of an iterative Newton-type method 
yields 

(7) Y„+i^k+\ = B{0ex,Yn;h„) +B{0i,„,Y„+i,k;hn) +A(F„+i,i+i - Y„+i,k) 

with some approximation to the Jacobian of B(<J>,m,I'n+i,jt;^n) and a predictor F„_|_i o. In 
the following we assume that (7) can be solved uniquely at least for sufficiently small 
To simplify the presentation, we assume further that all step sizes are constant, /i„ = h. 

For the approximation /;t there exist several common choices. If we choose Jj; to be the 
exact Jacobian d2B{<i>i,„.Y„^i i;;h), then we obtain the classical Newton iteration method 
for solving (5), with quadratic convergence. It will be denoted in the following as full 
Newton iteration. If we choose instead = d2B{^im,Yn,h), then we obtain the so called 
modified Newton iteration method, which is only linearly convergent. Here, Ji, is indepen- 
dent of the iteration number k, thus its computation is much cheaper and simpler than in 
the full Newton iteration case. The third and simplest possibility is to choose equal to 
zero. In this case we don't even have to solve a linear system for Yn+\.k+\. This iteration 
method is called simple iteration method or predictor corrector method. Its disadvantage 
is that for stiff systems it requires very small step sizes to converge. 

For most stiff problems and problems with additive noise, the use of semi-implicit meth- 
ods suffices. As will be demonstrated in Section 4 these have the advantage that less itera- 
tions are required to obtain the correct order of the underlying method. 



Here, 



^(i),/.„=^(Wi)-W(f„)=AW„ 




B(<?'«,i;;/J«) = -^0 + - ako(i'«) +/(i),/,„(l - /3)^i(i'„) 
+ (/(M),/>„-i3/fi)J {g\gx){Yn), 

B{(^im,Yn+\;hn) = hnag(){Yn+\) + I(\),h,fi g\{Yn+\) ■ 



^ex{»o)=hn{l-a), 
^e.(*l)=/(lM„(l-/3), 
«'e.([*l]l])=/(l,l),ft„-/3/('l),,„- 



<J»,m(»o) = ^«a, 

«'im(»l)=/(l),ft„/3 
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2. Stochastic B-series 

B-series, symbolized by B{^,xo;h), for SDEs were first constructed by Burrage and 
Burrage [2, 3] to study strong convergence in the Stratonovich case. In the following years, 
this approach has been further developed by several authors to study weak and strong con- 
vergence, for the Ito and the Stratonovich case, see e. g. [7] for an overview. A uniform 
theory for the construction of stochastic B-series has been presented in [7], in [8] this ap- 
proach has been used to construct order conditions for implicit Taylor methods. Following 
the exposition of these two papers, we define in this section stochastic B-series and present 
some preUminary results that will be used later. 

2.1. Some useful deiuiitions and preliminary results. 
Definition 2.1 (Trees). The set ofm + \-colored, rooted trees 

T = {0}U7bU7iU---U7;, 

is recursively defined as follows: 

a) : The graph = [0]/ with only one vertex of color I belongs to Ti. 

Let T = [ti . T2 , . . . . T^] / be the tree fanned by joining the subtrees Ti , T2 , . . . , T^r each by a 
single branch to a common root of color I. 

b) : IfxuX2,...,x^&Tthenx=[T:i,Vi,...,T^]i&Ti. 

Thus, 7/ is the set of trees with an /-colored root, and T is the union of these sets. 

Definition 2.2 (Elementary differentials). For a tree r gT the elementary differential is a 
mapping F{t) -.W^ defined recursively by 

a) : F(0)(xo)=xo, 

b) : F{»i){xo)=gi{xo), 

c) : //t = [ti , T2, . . . , Tk\i e 7] then 

F{T){xo)=gl''\xo){F{Ti){xo),F{T2){xo),...,F{T^){xo)). 
With these definitions in place, we can define the stochastic B-series: 
Definition 2.3 (B-series). A (stochastic) B-series is a formal series of the form 
B{(j),xo;h) = £ a(T) • (/»(T)(/i) • f (t)(xo), 

where ^ : T — >■ S := {{(p{h)}h>o ■ <Pih) : Q.^M. is Borel-measurable Mh > O} assigns to 
each tree a random variable, and a:T^Qis given by 

1 

a(0) = l, a(.,) = l, a(T = [Ti,...,T^],) = ^^ tYI^'^j)' 

'1 !'2! " ' rql 

where rj , r2, . . . , count equal trees among Ti , T2, • • • , Tk-. 

Note that a{x) is the inverse of the order of the automorphism group of T. 

To simplify the presentation, in the following we assume that aU elementary differentials 
exist and all considered B-series converge. Otherwise, one has to consider truncated B- 
series and discuss the remainder term. 

The next lemma is proved in [7]. 
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Lemma 2.1. IfY{h) = B{(j),xo;h) with </)(0) = 1 is some B-series and f e C°° (K^ , M'^"), 
then f{Y{h)) can be written as a formal series of the form 

(8) f{y{h)) = P{u)-n{u){h)-G{u){xo), 



ueUf 



where 



a) : Uf is a set of trees derived from T as follows: [0]/ € Uf, and if 

Ti,T2,...,Tk: e T\{9}then [Ti,T2,...,Tk:]/ e 17/, 

b) : Gmf)ixo) = fixo) and 

G(m = [Ti, . . . , T^]f){xo) = f^^^Xo) (F(Ti)(xo), . . . ,F{T^){xo)), 



c) : J3 ( [0]/) = 1 and J3 (m = [Ti , ... , tJ/) = — — 
w/jere ri,r2,...,rq count equal trees among Ti , T2, . . . , Tk:, 

d) : \l/<t>iMf) = 1 and \i/^{u=[T:u---,^K]f){h) =Uj=i(l>{^j){h). 

For notational convenience, in the following the /j-dependency of the weight functions 
and the XQ-dependency of the elementary differentials will be suppressed whenever this is 
unambiguous, so <I>(t)(A) will be written as <I>(t) and F(t)(xo) as F{t). 

The next step is to present the composition rule for B-series. In the deterministic case, 
this is e. g. given by [10], using ordered trees. The same rule applies for multicolored trees, 
as in the stochastic case. But it is also possible to present the result without relying on 
ordered trees, as done in [8], and this is the approach that will be used in the following. 

Consider triples (t, i?, (b) consisting of some T € T, a subtree i5 sharing the root with T, 
and a remainder multiset (O of trees left over when j> is removed from t. We also include 
the empty tree as a possible subtree, in which case the triple becomes (t,0, t). 





Example 2.1. Two examples of such triples are 

\ 

and 

y V ) 

So, for the same T and i> there might be different CD's. 

We next define ST{t) as the set of all possible subtrees of T together with the corre- 
sponding remainder multiset ft), that is, for each T G r\0 we have 

5r(.,) = {(0,.,),(.,,0)}, 

5r(T=[Ti,...,T^]/) = : I? = [t?i,...,t?^];, (0 = {(0i,...,(0k}, 

coi) e STizi), J = 1, . . . , If } u (0,t). 

We also have to take care of possible equal terms in the formula presented below. This is 
done as follows: For a given triple (t, t>, ffl) write first ^5 = [tJi , . . . , = , . . . , 4'],, 
where the latter only expresses that is composed by q different nonempty trees, each 
appearing r,- times, hence Y!i=i = . Let T = [ti , . . . , Tk-];. For i = l,...,q, each t?,- is a 
subtree of some of the T, 's, with corresponding remainder multisets coj. Assume that there 
are exactly different such triples (t/a,!?,-,©^) each appearing exactly r,-,^ times so that 
Lfii = ^i- Finally, let S^G (Ohe the distinct trees with multipUcity Sk,k = l,...,po,of 
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the remainder multiset which are directly connected to the root of T. Then, t can be written 
as 

l^-' i- — L''l ) • • • > "po ) ''11 ) • • • > ''Ipi '■ql 5 • • • > J( ~ L'-l ) • • • ) ''2 \h 

where the rightmost expression above indicates that T is composed by Q different trees 
each appearing times. 

With these definitions, we can state the following theorem, proved in [8|: 

Theorem 2.1 (Composition of B-series). Let (l>x,(l>y ■ T and <^>«(0) = 1. Then the 
B-series B{(j)x,xo;h) inserted into s/z) is again a B-series, 

B{(j)y,B{(j)x,xo;h);h) ^ B{(j)xO (j)y,xo;h), 

where 

(l>,(B)6Sr(T) \ 5€C» / 



7(0,0,0) = 1, and 

Ril - Rnl AA . K 



Sl'. Spg.ru- f^qpq- i=lk=l 

for T given by (9). 

The combinatorial term 7 gives the number of equal terms that will appear if the com- 
position rule using ordered trees is preferred. 

In general, the composition law is not Unear, neither is it associative. It is, however, 
linear in its second operand. Further, if both <j>x{9) = ^> (0) = 1> then the composition law 
can be turned into a group operation (Butcher group, see [6, 10, 11] for the deterministic 
case): The inverse element <^>"' (t) can be recursively computed by 



(10) {(l)o(j>-'){T)=e{T) 



1 ifT = 0, 
otherwise, 



and associativity is proved by 

(11) B{<j)^,B{<j)y,B{<j)x,xo;hy,hy,h) =B{<j)xo{<j)yO<j)^),xo;h) =B{{<j)xO<j)y)o<j)^),xo;h). 

This holds even if 0^(0) ^ 1. 

The next result will be needed for the investigation of modified Newton iterations. 

Lemma 2.2. If(j)x(V)) =Owe have 

d2B{<l>y,xo;h)B{(px,xo;h) =B{<(>x*<l>y,xo;h), 

where the bi-Unear operator * is given by 

[0 //T = 0, 

(12) {^^^^y){x)=l J. y{x,-&,{5})-<i>y{-&)<i>x{d) otherwise, 
with 

SP{t) = {{■&,(o) G ST{t) : (0 contains exactly one element 8}. 



8 



KRISTIAN DEBRABANT AND ANNE KVjERN0 



Proof. Written in full, the statement of the theorem claims that 
i^eTSeT 

(13) = £ a(T)( £ r{^,^,{5})-(^y{&)U5)]-FiT). 

T€r\{0} V('?,{5})eSP(T) / 

This is true if 

(14) {dF{^)F{5))= £ j3(T,t>,5)F(T) 
and 

(15) a(t?)a(5)/3(T,t?,5) = a(T)7(T,i?,{5}), 

where 5) is the set of all t's constructed by attaching 5 to one of the vertices of i?. 

We will prove this by induction. 

First, let = 0. Since dFiQl)F{5) = F{5) we have T = 5 and (14) and (15) are trivially 
true with j3(T,0,T) = 1. Now, let i} = Then dF{i})F{8) = g\F{5) = F{[S\i). As 
= {[5]/} this gives j3([5]/,»/,5) — 1, and again (15) is trivially true. Finally, let 

^ = W,\...,-&l',...,¥,'']i 

with distinct trees t?i , . . . , t?^. Then 

ri times times 

aF(^?)F(5)=g|'^*+')(F(5),F(di),...,F(di),...,F(^,),...,F(d,)) 

r, — 1 times 

+ 1 n^l'"' • • .,dF{^i)F{8)',F{^i)S..,F{-&ii . . . ,F{d,)) , 

i=\ 

where fc^ = Y^^i rt, so T e ■!a?'(t?, 5) is either 

T = [5,^?[^...,^5^^...,t?;»], with J3(T,1?,5) = 1 

or 

T = [^[' ,...,Ti, with Ti e 5) and j3 (T, tJ, 5) = r,j3 (t,-, tJ,-, 5). 

In the first case, if 5 = t?; for some j, then a(T) = a(tJ)a(5)/M and /(t, tJ, {5}) = M 
with M = ■ + 1. Otherwise the same is valid with M = I. So (15) holds. In the second 
case, assume that our induction hypothesis (15) is true for all T, e (i?,, 5). We obtain 

«(T) = ^^aW and r(T,t?,{5})=Mr(T,-,d,-,{5}) 

with M = ■ + 1 if T, = for some j and M = 1 otherwise. It follows that 

a(T)r(T,t?,{5}) = r,|M7(T,-,i?,-,{5})a(t?) = a(t?)a(5)r,j3(T,-,d,-,5) 
= a(i?)a(5)j3(T,t?,5). 

□ 
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2.2. B-series of the exact and the numerical solutions. From the results of the previous 
subsection, it is possible to find the B-series of the exact and numerical solutions. Here, 
the proofs are only sketched, for details consult [7, 8]. 

Theorem 2.2. The solution X {to + h) of (4) can be written as a B-series B{(p,xo;h) with 

(p(0) = l, (p{»i){h)=Wi{h), (p{T=[ti,...,T^]i){h)= [''ll(p{Tj){s)^dWi{s). 

JO j=i 

Proof. Write the exact solution as some B-series X{tQ + h) = B{(p,xo;h). As 9(0) s 1, 
apply Lemma 2.1 to gi(X{to + h)) and Definitions 2.1, 2.2, and 2.3 to obtain 

(16) giiB{cp,xo;h)) = £ a(T)-(?>;(T)(/i)-F(T)(xo) 

TETi 

in which 

'1 ifT = »/, 

<P/'(T) W = ^(^.^(^^ if T = [Ti, . . . , T,]i G Ti. 

Insert this into the SDE (3) and compare term by term. □ 
Theorem 2.3. The numerical solution Yi given by (5) can be written as a B-series 

Yi=B{^,xo;h) 

with <J> recursively defined by 

(17a) <I>(0) = 1, 

(17b) «I>(t) = 0ex{t) + {^O0i^){T). 

Proof Write = B{(i>.XQ;h) and insert this into (5). As 3>(0) = <^exi®) + '^im{®) = 1, 
apply Theorem 2.1, and compare term by term. □ 

To study the consistency of the numerical methods, we need to assign to each tree an 
order: 

Definition 2.4 (Tree order). The order of a tree x gT respectively uGUf is defined by 

K 

p(0) = 0, p(m = [Ti, . . . , tJ/) = £p(tO, 

(=1 

and 

IX /X fl for 1^0, 

p(T=[Ti,...,Tj/) = 2^p(Ti) + < 1 ■ , 

" I 2 otherwise. 

In Table 1 some trees and the corresponding values for the functions p, a, and (p are 
presented. 

To decide the weak order we will also need the B-series of the function /, evaluated 
at the exact and the numerical solution. From Theorems 2.2 and 2.3 and Lemma 2.1 we 
obtain 

f{X{to + h)) = £ J3(m) • W<p{um ■ G{u){xo), 

ueUf 

m) = £ • ¥^{u)ih) ■ G{u){xo), 
ueUf 
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P(t) a(T) 



.h if/ = 



1' 



[^(2,0,1) W 



2 i /oViW^^J. ^ J 2^1,1,0) (S) 

(2/(1,1,0) +%o) (1) 



5 Jo do' w^2(^2)'*^;wi(^2)) 





[4^(2,2,1,1,0) +2/(2,1,2,1,0) 


+ 2/(1,2,2,1,0) 


(S) 




1 4/(2,2,1,1,0) + 2/(2,1,2,1,0)' 


+"2/(1,2,2,1,0) 






[ +2/(0,1,1,0) + 2/(2,2,0,0) + 


'^(1,0,1,0) +-^(0,0,0) 


(I) 



Table 1 . Examples of trees and corresponding functions p (t), 05(t), 
and (P(t). The integrals (p(T) are also expressed in terms of multiple 
integrals 7( ) for the Stratonovich (S) and /( ) for the Ito (I) cases, see 
[15] for their definition. In bracket notation, the trees will be written as 
[[•2\o\i, [•i,»i]o, and [•!, [#2, •2] i]o, respectively. 



with 

K 

MMf) = h v^,p(M = [ti, . . .,T^]f) = n (Pi^j), 

and 

K 

V/<i.([0]/-) = 1, W<s>{u = [ti,...,tJ^) = n^('^i)- 

;=i 

One can show [15, 5, 9] that E\i/^{u){h) = ^(F^")) Vm e Uf and 9(t)(/i) = 
Vt e T, respectively, where especially in the latter case the ^(•)-notation refers to the 
L2(n)-norm and h^O. 

In the following we assume that also method (5) is consistent with the definition of the 
tree order, i.e. that it is constructed as usual such that E\j/^{u)(h) = &{hP^"^)\/u G Uf and 
3>(t)(/i) = ^(/jP^"^)) Vt e T, respectively. These conditions are fulfilled if for t € T and 
e N = {0, 1, . . . } it holds that {^ex{r)f = and (3>,>„(t))^* = &{h^'P^^'^). 

3. B-SERIES OF THE ITERATED SOLUTION AND GROWTH FUNCTIONS 

In this section we will discuss how the iterated solution defined in (7) can be written in 
terms of B-series, that is 

Assume that the predictor can be written as a B-series, 
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1 ^ V 

{)(T)=3,t(T)=5(T) = l; F,(T)=4,r(T) = 3,J)(T)=2; f,(T) = r(T) = = 3 

Figure 3 . Examples of trees and their growth functions for simple (f)), 
modified Newton (r) and full Newton (0) iterations. 

satisfying <t>o(0) = 1 and 00(1) = ffih^^'^^) Vt G T . The most common situation is the use 
of the trivial predictor Y\ = xq, for which <t>o(0) = 1 and <t>o(T) = otherwise. 

We are now ready to study each of the iteration schemes, which differ only in the choice 
of Jk in (7). In each case, we will first find the recurrence formula for ^^(t). From this we 
define a growth function ^(t): 

Deftnition 3.1 (Growth function). A growth function g : T — >■ N w a function satisfying 

Oj^(t) = <I)(t) Vt e r with q{t) < k 
'■^^^ => <&;t+i (t) = <I>(t) yTeTwithQ{T)<k+l, 

forallk>0. 

This result should be sharp in the sense that in general there exists T ^ with «&o(t) ^ 
<J>(t) and ^k{^) ^ 4>(t) when k < g{T). From Lemma 2.1 we also have 

fiYi,k)^ £ j3(«)-VA<i.,(m)-G(m)(xo) 
ueUf 

with 

K 

V^Mf) = 1' V^*,(m = [Ti,...,T^]/) = Yl^kiTj), 
where J3(m) and G{u){xo) are given in Lenmia 2.1. This implies 

(19) <^<i.t(M) = r<i.(T) VM=[Ti,...,T^]/et//with0'(M)=rnax0(T,)<L 

As we will see, the growth functions give a precise description of the development of 
the iterations. However, to get applicable results we will at the end need the relation 
between the growth functions and the order. These aspects are discussed in the next section. 
Examples of trees and the values of the growth functions for the three iteration schemes 
are given in Figure 3. 

3.1. The simple iteration. Simple iterations are described by (7) with = 0, that is 

(20) Y„+^,k+^=Bi<^>e.,Yn■,h)+B{<^>i,„,Y,+ ^y,h). 

By Theorem 2.1 we easily get the following lemma, where, as in the following, all results 
are valid for all / = 0, . . . ,m: 

Lemma 3.1. IfYi o = B(4>o,xo;/i) then Y^^^ = B{<t>k,XQ;h), where 

<I'i+i(0)=l, 

^k+l {t:) = ^ex(T) + {^k ° ^<m) 
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The corresponding growth function is given by 

()(0)=O, ()(.,) = 1, f)([Ti,...,T^];) = l+max()(Tj). 

The function f) (t) is the height of T, that is the maximum number of nodes along one 
branch. 

3.2. The modified Newton iteration. In this subsection we consider the modified Newton 
iteration 

Yn+l,k+l=B{^ex,Yn,h)+B{^i„^,Y„+i^k;h) 
^^^^ + d2B{^i,n,Yn\h){Yn+l,k+l - Yn+l,k). 

The B-series for Yi^k and the corresponding growth function can now be described by the 
following lemma: 

Lemma 3.2. IfYx o = B{<^o,xo;h) then Yij^ = B{<^k,xo;h) with 
The corresponding growth function is given by 

r(Ti) if K=l, 



t(0)=O, r(.,) = l, r(T=[Ti,...,Tj;): 



H-maxr(Ty) if k>2. 



The function r(T) is one plus the maximum number of ramifications along any branch 
of the tree. 

Proof. Theorem 2.1 and Lemma 2.2 imply (22). We next prove that r is the appropriate 
growth function. If r(T) = then T = and 4>o(t) = <I>(t). Assume now that <I>a;(t) = <I>(t) 
Vt with r(T) < k. If G ST{t)\SP(t) with r(T) <k+l, then V5 e © it holds 

r(5) < r(T) and therefore 3>fc(5) = 3)(5). So, by (22) we have Vt with r(T) <k+l 

3'*+l(T)=3'e.(T)+ 7(T,1?,©)-<f»(t?)n'f(5) 

(j>,(B)Gsr(T)\s/'(T) dea 
+ £ r(T,i?,{5})-4>™(t?)4>,+i(5), 

and by induction on the number of nodes of these trees we obtain that ^k+i ('t) = ^{'^) 
withr(T) <^+l. □ 

3.3. The full Newton iteration. In this subsection we consider the full Newton iteration 
(7) with Jk = d2B{'i>im,Yi^k'->h). Extending the *-operator to the case when its first operand 
does not vanish on the empty tree by 

{(l>x * = {{(l>x - 0^(0)e) * + 0^(0)«|>y(T), 

it follows that the B-series for Yi^k and the corresponding growth function satisfy: 
Lemma 3.3. IfYi^ = B((^o,xo;h) then Yi^k+i = B{(^k+i,xo',h) with 

<^k+im=h 

(23) 
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The corresponding growth function is given by 

t)(0)=O, £)(•,) = 1, 

r IN fmax'^L, ()(t;) if 7=1, 

[max}Lit)(T;) + l if 7>2, 

w/iere y the number of trees in Ti,...,Tk satisfying t)(T,) = maxj^j J)(Ty). 

The function 5 is called the doubling index of T. 

Proof. Writing the iterations in terms of B-series, we get 

(24) B(3);fc+i,xo;ft)=B(4>e„xo;/i)+B(4>,„„yi,ft;ft) + a2B(<I>,„„I'i,fc;/j)B(A3>;fe,xo;/j) 
with A^i.(T) = 3>fc+i(T) -3>fc(T). I^txo=B{%\Yi^},;h) so that 

The use of Lemma 2.2 followed by the use of Theorem 2. 1 give the following result: 
d2B{^i,„,Y^^k\h)B{^l^ oA^k,yi,k\h)=B{{^-' oA^k) *^i,n,y\.k-^h) 

- B(4>i o ( (4>^ 1 o A^k) * ^im),xo; h) . 
The operator * is bilinear and o is linear from the right, thus 

and the first part of the theorem is proven by (24). We will now prove the second part. 
Assume that 4><;(t) = (i>k+i (t) = «I>(t) for all T satisfying 5(t) < k. This is true for = 
and T = 0. Let ^'|( = <I>^' °^k+i and notice that by the assumption above, ^/:(t) equals 
the unit element e(T) if 0(t) < k. Consider a tree T where c)(t) = k+1. For this tree we 
obtain 

(25) + £ r(T,t?,aj)4>,,„(t>)ni'-t(5) = Cl'i*<I>,m)(T) 

(i?,(B)€sr(T)\si'(T) 8e<o 

since the last sum of (25) disappears: For each (t3-,co) G ST{t) \ SP{t) (if any) there is at 
least one 5 G ft) satisfying d{5) < k and thereby ^k{S) = 0. In this case we obtain 

(4>iO ((4>-i o<I>^+i) *<I>,„,))(t) = (OioO^ 1 oO^^+i o <!>,„,) (t) = {^k+io^i^){T) 

by (11), so that 

The theorem is completed by induction on the number of nodes of T and on k. □ 
4. General convergence results for iterated methods 

Now we will relate the results of the previous section to the order of the overall scheme. 
We have weak consistency of order p if and only if 

(26) EY^{u){h)=EY^{u){h) + ^{hP+'^) W GUf with p{u) < p + ^ 

((26) slightly weakens conditions given in [16]), and mean square global order p if [4] 
4>(t)(/j) = (P(t)(/j) + ^(/i^+2) VTe7'withp(T)<p, 
E^{T){h) = E(p{T){h) + ^{hP+^) VTeTwithp(T)</7 + i, 
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(27) Eva<j,^(m) =Eva^(m) Vm e f// with p (m) < + - 



and all elementary differentials F(t) fulfill a linear growth condition. Instead of the last 
requirement it is also enough to claim that there exists a constant C such that < 
C V>' e M™, 7=0,... ,M, and all necessary partial derivatives exist [3]. 
Then, the order of the iterated solution after k iterations is if 

1 
2 

in the weak convergence case respectively 

^k{T:) = <P(t) Vt e r with p(t) < qk, 

(28) 1 
Ea>fc(T)=E(p(T) VTe7'withp(T)=gft+- 

in the mean square convergence case. 

In the following, we assume that the predictors satisfy the condition 

(29) <I>o(t) =4>(t) Vtg r with0(T) 

where % is chosen as large as possible. In particular, the trivial predictor satisfies % = 0. 
It follows from (18) and (19) that 

(30) a>ft(T)=3>(T) VTeTwith0(T) <%+A;, 
as well as 

(31) V**!") = V*(") € Uf with 0'(m) <% + k. 

The next step is to establish the relation between the order and the growth function of a 
tree. We have chosen to do so by a maximum growth function, given by 

(32) ^(^) =m^{fl(T) : p(t) < ^} = f^{^{u) ■.p{u)<q}. 

With this definition, by (31) respectively (30), the conditions (27) respectively (28) are 
fulfilled for all u of order p{u) < vtan{q]^,p) respectively aU t of order p{x) < min(^^,/7) 
if 

(33) '^{qk + \)<% + k. 

Let CT and C Uf be the set of trees with an even number of each kind of stochastic 
nodes. E. g. from [9] we have 

e^(t)=o if T^r^, 

(34) 

E\i/^{u)=0 if u^Uf. 

Thus, if the method is as usual constructed such that also Vm,n € N and Vti^, €T, i = 
1,...,OT, Vt2j e 7", 7 = l,...,n. 



(35) E YlYl^ex{Tl,i)^im{T2j) ]=0 if £p(Ti,,) + £ P{T2j) i N, 

V<=17=1 / i=l 7=1 

then in (33) qk-\- \ can be replaced by \_qi + \\ . 

The results can then be summarized in the following theorem: 

Theorem 4.1. If (35) is fulfilled, then the iterated method is of weak respectively mean 
square order q^ < p after ^{[qt + j\) iterations, otherwise after ^{qi^ + j) — % 
iterations. 
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1 V 

^(,,3 ^t,3 ^0,3 

Figure 4. Minimal order trees with 0(t) = 3. The sets =^,3 consist of 
all such trees with only stochastic nodes. 

Our next aim is to give explicit formulas for the maximum growth function. Let us start 
with the following lemma. 

Lemma 4.1. For k>l, 

^r)=k p(T)>2*-i-i 
The same result is valid for i)'{u), x' (u), and Z)'(m). 

Proof. Let ^f,^^^, 5^^^^, and ^-a^k be sets of trees of minimal order satisfying f)(T) = A: Vt e 
=^f).*' '-('^) = kVx G ,%,k, and X){x) = kVT E (see Figure 4), and denote this minimal 
order by P(,^i, p^^k^ and pxi^k- Minimal order trees are build up only by stochastic nodes. It 
follows inomediaitely that ^f,,i = = = {•/ : / > 1}. Since p(»/) = 1/2 for / > 1, 
the results are proved for A: = 1 . It is easy to show by induction on k that 

1 k 

^^,k = {['^]i- i^e /> 1}, p[,_fe = p(,,^_i + - = 

(36) ^r,/t = {[•/,, : TG /l,/2 > 1}, Pa=Pr,*:-l + l=^-^, 

=^0,* = {[ti,T2]/ : Ti,T2 e =3^£,,i-i, / > 1}, pi,,k = 2po,t-i + ^ = 2*-i - i. 

For each g being either f), r, or 0, the minimal order trees satisfying 0'(Mg,ft) = A: are Mg ,t = 
[Tg i]/ with Tg jt e =3g.jt, which are of order p(Tg^yt). □ 

Now we can prove the following corollary. 

Corollary 4.1. For q> jwe have 

{2q for simple iterations, 

[q + j\ for modified Newton iterations, 

[log2 + 5 ) J +1 for full Newton iterations. 

Proof. The minimal order trees are also the maximum height / ramification number / dou- 
bling index trees, in the sense that as long as p(Tg_jt) <q< p{Tg^k+i) there are no trees of 
order q for which the growth function can exceed k. □ 

For some methods, these results can be refined. We call a method semi-implicit, if 
(t) = Vt ^ Jo (remember that Tq is the set of trees with a deterministic root). Then, 
by Lemmas 3.1, 3.2, and 3.3 we obtain the following lemma: 
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p 


simple iter. 


mod. iter. 


full iter. 


1 

2 


2(1) 


1 


1 


1 


2(1) 


1 


1 


1^ 


4(2) 


2 


2 


2 


4(2) 


2 


2 




6(3) 


3(2) 


2 


3 


6(3) 


3(2) 


2 



Table 2. Number of iterations needed to achieve order p when us- 
ing the simple, modified or full Newton iteration scheme in the Ito and 
Stratonovich case for strong or weak approximation, provided (35) is 
fulfilled. In parentheses, the numbers for semi-implicit methods are 
given. 



Lemma 4.2. For semi-implicit methods, the corresponding growth functions are given by 



1 if I > 0, 



ts(0)=O, Xs{»l) = l, rs(T= [Ti,...,tJ; 



ts(Tl) 



if l>0, 

if / = 0,ff=l, 



1 + max ts (t,) if / = 0, fc > 2, 



0^(0) = 0, c)s(«/) = l, Z)s(t=[ti,...,t^],) = < 



1 



«/ />o, 
if / = 0,7=1, 



niaxc)5(T;) 
max05(T^) + l !/ l = 0,Y>2, 



w/jere 7 is the number of trees in Ti,...,Tk satisfying (Tj) = maxj^j ("^y )• 

This implies immediately: 
Lemma 4.3. For k>l, 

> p(T)>A;-i 



rs(T)=A; 

0<r(T) =fc 

r/ie iame reiM/? is valid for i)s'{u), ts'{u), and ds'{u). 
Corollary 4.2. For semi-implicit methods we have for q> j 



P{r)>-k-l, 
p(T)>^2^-l. 



[q + j] for simple iterations, 

[| + 1 )J for modified Newton iterations, 

[log2 ^Y^-J + 2 for full Newton iterations. 



For the trivial predictor. Table 2 gives the number of iterations needed to achieve a 
certain order of convergence, both in the general and in the semi-implicit case. 
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For the sake of completeness, we also give the corresponding results for (deterministic) 
Taylor methods applied to deterministic problems. Note that in this case, (35) is automati- 
cally fulfilled. 

Lemma 4.4. Suppose that the considered problem is purely deterministic, i. e. m = in 
(3). Then, for k> 1, 

^ P{'c)>k, 

x{T)=k ^ p(t)>2A;-1, 

J)(t)=A: ^ p(t)> 2*^-1. 

The same result is valid for [)'(«), x' (u), and 9s' (m). 

Corollary 4.3. For deterministic problems, we have for q&^, q>\ 

q for simple iterations, 

^(^) = S L^T^J for modified Newton iterations, 

, [logiC^ + 1)J for full Newton iterations. 

5. Numerical examples 

In the following, we analyze numerically the order of convergence of several stochastic 
Taylor methods in dependence on the kind and number of iterations. 

As first examples, we apply the semi-implicit Milstein method [15], denoted by SIM 
and given by (6), the implicit Milstein-Taylor method [17], denoted by IM and given by 

Yn+l = Y„ + hgQ{Y„+i)+I(i)gi{Yn+l)- (7(1,1) +/!)k'i5i](y„+i), 

both of strong order 1.0, and the semi-impUcit strong order 1.5 Taylor method due to 
Kloeden and Platen [15, 17], denoted by SIKP and given by 

F„+i =i'„+/igo(l'«+i)+/(i)^i(l'«)+7(ia)k'igi](i;)-7(o,i)kogi](i'«) 
- \h^[8'oSo + \gU\\ +7(0,1) [^'i^o + \g'[8\] {Yn) 
+ h\,\Aghi+g'ig'i\{Yn), 
to the non-Unear SDE [15] 

(37) dX{t) = (^X{t) + ^Jx{tY + \^ dt + ^Jx{tY + \dW{t), X(0) = 0, 

on the time interval / = [0, 1] with the solution = sinh(f + W{t)). With each method, 
the solution is approximated with step sizes 2~^\...,2~^^ and the sample average of 
M = 4000 independent simulated realisations of the absolute error is calculated in order to 
estimate the expectation. 

The results at time ? = 1 are presented in Figure 5, where the orders of convergence 
correspond to the slope of the regression lines. As predicted by Table 2 we observe strong 
order 1 .0 for one simple or one (modified) Newton iteration of the semi-implicit Milstein 
method; and no convergence for one simple iteration but strong order 1 .0 for two simple or 
one (modified) Newton iteration of the implicit Milstein-Taylor method. The semi-implicit 
strong order 1 .5 Taylor method yields strong order 1 .0 for one and strong order 1 .5 for two 
simple or modified Newton iterations. 
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-SIM(SIl): p=0.99 
-SIM(MIl): p=1.00 
IM(SIl): p=-0.01 
-IM(SI2): p=1.00 
-IM(MIl): p=0.98 












^SIKP(SIl): p=0.99 
-o- SIKP(SI2): p=1.53 
^SIKP(MIl): p=1.00 
- B -SIKP(MI2): p=1.50 



-13 
log, h 



-13 
log, h 



(a) Semi-implicit (SIM) and implicit (IM) Milstein (b) Semi-implicit strong order 1.5 Taylor (SIKP) 
method method (the results for two simple or modified New- 

ton iterations nearly coincide) 

Figure 5. Error of several Taylor methods applied to (37) with up to 
two simple (SI) and modified Newton (MI) iterations 



Next, we apply the semi-implicit weak order two Taylor scheme due to Platen [15], 
denoted by SIW and given by 



y„+i =Y„ + hgQ{Y„+i) + 7(1)^1 (y„) +/(i,i) fe'i^i] (i;) 



1 



to SDE (37). Here, we choose as functional /(jc) = p(arsinh(jic)), where p{z) — 6z^ + 8z 
is a polynomial. Then the expectation of the solution can be calculated as 

(38) E(/(X(f)))=f'-3f2 + 2f . 

The solution E{f{X{t))) is approximated with step sizes 2 2-^ and M = 4 • 10*^ 

simulations are performed in order to determine the systematic error of SIW at time t ~ I. 
The results with one or two simple or modified Newton iteration steps are presented in 
Figure 6. According to Table 2 we expect approximation order one for one iteration and 
order two for two iterations, which is approved by Figure 6. Finally, we apply the fully 
implicit strong order 1.5 Taylor scheme given in [8], 

1 1 K ^ , 1 2 , 

Y„+l =Y„ + -/(1)§1,„+1 + -hgO,n+l + 2 (^(1,1) +'')^L«+l^l-n+l + 4'^ gO.n+lgO.n+l 

+ ^h^gO,n+l{gi,»+i,gUn+l) + ^I{l)gl + ^hgQ - + g'lgl 
+ 1^ (^0,1) -'^(l,0))/l§0- t(Vi) -'f(KO)ko^l 



1 7 

+ 1 2^(0,1) -4^(1) -2/, 



(1,1,1)^ ^ifei'^i)- 



/!/(l)+2/(i 1 1) i^i^i^i 



h^g'ogo - h^g"{go^gi)~ ^f^^gigogi 



^Vi/i^ 



-h^a&x.gi) 



h^g'igWigi - ^h^g'ig'l(gugi) - lh^g'({g'ig,gi)- ^h^g'iigugugi, 
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-«-SlW(SIl): p=0.90 




^S1W(SI2): p=1.81 




-o- SIW(MU): p=0.90 




^SIW(MI2): p=l.81 



log^h 



Figure 6. Error of the weak second order (semi-) implicit Platen 
method applied to (37) with one or two simple (SI) or modified Newton 
(MI) iterations (the results for one respectively two simple and modified 
Newton iterations coincide) 



(here we used the abbreviations g/„+i — gi{Y„+i) and gi — gi{Y„)), which is denoted by 
FIT, and the semi-implicit strong order 1 .5 scheme SIKP to the system of non-linear SDEs 

(39) 

dXi{t)= Qxi (f ) + ^JXi {tf +X2{tf + 1^ t/f+(sin(Xi(f))+2sin(X2(f))) t/W(f), 

dX2{t)^ Q^i(r) + \/^2(r)^ + l) t/f+(cos(Xi(f)) + 3cos(X2(f))) t/W(f), 

Xi(0)=0, X2(0)=0, 

again on the time interval / = [0, 1]. The solution is approximated with step sizes 2^" , ... ,2 
and the sample average of M = 4000 independent simulated realisations of the absolute er- 
ror is calculated in order to estimate the expectation. As here we do not know the exact 
solution, to approximate it we use SIKP with two simple iterations and a step size ten times 
smaller than the actual step size. 

The numerical results at f = 1 are presented in Figure 7. Again, the orders expected 
according to Table 2 are confirmed. 

6. Conclusion 

For stochastic implicit Taylor methods that use an iterative scheme to approximate the 
solution, we derived stochastic B-series and corresponding growth functions. From these, 
we deduced convergence results based on the order of the underlying Taylor method, the 
choice of the iteration method, the predictor, and the number of iterations, for Ito and 
Stratonovich SDEs, and for weak as well as strong convergence. The convergence results 
are confirmed by numerical experiments. From a practical point of view, this theory might 
lead to the construction of more efficient numerical schemes for SDEs. But we also like 
to point out that the similarities of the iteration dependent growth functions g for a range 
of problems (ODEs, DAEs, and SDEs) and underlying methods (Runge-Kutta methods, 
implicit Taylor methods) indicate an underlying structure that could well be investigated in 
a more general fashion. In spite of this, the number of iterations needed to obtain the order 
of the underlying implicit Taylor method does not depend on whether the SDE is of Ito or 



20 



KRISTIAN DEBRABANT AND ANNE KV^RN0 



-2 
-2.5 



-3 




-15 -14 -13 -12 -II 



log h 

(a) Semi-implicit strong order 1.5 Taylor (SIKP) 
method (the results for two simple or modified New- 
ton iterations coincide) 







FIT(SIl): p=-0.01 












-<-FIT(SI2): p=1.01 













^FIT(SI3): p=1.12 












-0- FIT(SI4): p=1.49 
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-«-FIT(MI2):p=1.50 
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(b) Full implicit strong order 1.5 Taylor (FIT) method 
(the results for four simple or two modified Newton 
iterations nearly coincide) 



Figure 7. EiTor of SIKP and FIT applied to (39) with different num- 
bers of simple (SI) and modified Newton (MI) iterations 

Stratonovich type. This is in contrast to the results obtained for Runge-Kutta methods for 
SDEs, for which usually less iterations are needed in the Stratonovich case [7]. The reason 
for this is that in the latter case certain error terms have vanishing expectation even if they 
do not vanish themselves. This is not the situation for implicit Taylor methods. 
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